/*
 * Matrix.cpp
 *
 *  Created on: Dec 24, 2014
 *      Author: Phong
 */

#include "Matrix.h"
#include <cblas.h>
#include "Utils.h"

Matrix::~Matrix() {
    delete[] data;
}

Matrix::Matrix(int rows, int cols) {
	data = new real[rows*cols];
	this->rows = rows;
	this->cols = cols;
	length = rows * cols;
}

// vector
Matrix::Matrix(int rows) {
	data = new real[rows];
	this->rows = rows;
	this->cols = 1;
	length = rows;
}

Matrix::Matrix(int rows, int cols, real* data) {
	this->rows = rows;
	this->cols = cols;
	this->length = rows * cols;
    this->data = new real[length];
    memcpy(this->data, data, sizeof(real)*length);
}

Matrix::Matrix(int rows, real* data) {
	this->rows = rows;
	this->cols = 1;
	this->length = rows;
    this->data = new real[length];
    memcpy(this->data, data, sizeof(real)*length);
}

void Matrix::save(ofstream &f) {
    f << rows << " " << cols << endl;
    for (int i = 0; i < length; i++)
        f << data[i] << " ";
    f << endl;
}

Matrix* Matrix::load(ifstream &f) {
    int rs, cs;
    f >> rs >> cs;
    Matrix* x = new Matrix(rs, cs);
    for (int i = 0; i < x->length; i++)
        f >> x->data[i];
    return x;
}

void Matrix::print() {
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++)
            cout << data[i*cols + j] << " , ";
        cout << endl;
    }
    cout << endl;
}

real Matrix::get(int r, int c) {
	return data[r*cols + c];
}

void Matrix::put(int r, int c, real value) {
	data[r*cols + c] = value;
}

Matrix* Matrix::dup() {
	Matrix* ret = new Matrix(rows, cols);
	memcpy(ret->data, data, sizeof(real)*length);
	return ret;
}

void Matrix::copy(Matrix* other) {
	delete[] data;
	rows = other->rows;
	cols = other->cols;
	length = rows * cols;
	data = new real[length];
	memcpy(data, other->data, sizeof(real)*length);
}

void Matrix::errorIfNotSameSize(Matrix* other) {
	if (rows != other->rows || cols != other->cols)
		Utils::error("size not match");
}

void Matrix::fill(real value) {
	for (int i = 0; i < length; i++)
		data[i] = value;
}

// y <- alpha * x + y
Matrix* Matrix::addi(real alpha, Matrix* x) {
	errorIfNotSameSize(x);
	for (int i = 0; i < length; i++) {
		data[i] += x->data[i] * alpha;
	}
	return this;
}

// y <- x + y
Matrix* Matrix::addi(Matrix* x) {
	return this->addi(1, x);
}


// y[i] <- y[i] + a
Matrix* Matrix::addi(real a) {
	for (int i = 0; i < length; i++) {
		data[i] += a;
	}
	return this;
}

// z <- alpha * x + y
Matrix* Matrix::add(real alpha, Matrix* x) {
	errorIfNotSameSize(x);
	Matrix* z = new Matrix(rows, cols);
	for (int i = 0; i < length; i++) {
		z->data[i] = data[i] + x->data[i] * alpha;
	}
	return z;
}

// y <- x + y
Matrix* Matrix::add(Matrix* x) {
	return this->add(1, x);
}


// z[i] <- y[i] + a
Matrix* Matrix::add(real a) {
	Matrix* z = new Matrix(rows, cols);
	for (int i = 0; i < length; i++) {
		z->data[i] = data[i] + a;
	}
	return z;
}

// z[i] <- x[i] * y[i]
Matrix* Matrix::mul(Matrix* x) {
    errorIfNotSameSize(x);
    Matrix* z = new Matrix(rows, cols);
    for (int i = 0; i < length; i++) {
        z->data[i] = data[i] * x->data[i];
    }
    return z;
}


// y[i] <- x[i] * y[i]
Matrix* Matrix::muli(Matrix* x) {
	errorIfNotSameSize(x);
	for (int i = 0; i < length; i++) {
		data[i] *= x->data[i];
	}
	return this;
}

// y[i] <- a * y[i]
Matrix* Matrix::muli(real a) {
	for (int i = 0; i < length; i++) {
		data[i] *= a;
	}
	return this;
}

// z[i] <- a * y[i]
Matrix* Matrix::mul(real a) {
	Matrix* z = new Matrix(rows, cols);
	for (int i = 0; i < length; i++) {
		z->data[i] = data[i] * a;
	}
	return z;
}

Matrix* Matrix::divi(real a) {
	return this->muli(1/a);
}

Matrix* Matrix::divi(Matrix* x) {
    errorIfNotSameSize(x);
    for (int i = 0; i < length; i++)
        data[i] /= x->data[i];
    return this;
}

Matrix* Matrix::div(real a) {
	return this->mul(1/a);
}

Matrix* Matrix::sqrt() {
    Matrix* z = new Matrix(rows, cols);
    for (int i = 0; i < z->length; i++)
#ifdef __DOUBLE__
        z->data[i] = sqrtl(data[i]);
#else
        z->data[i] = sqrtf(data[i]);
#endif
    return z;
}


Matrix* Matrix::getColumn(int col) {
	if (col >= cols) Utils::error("exceed num of cols");

	Matrix* ret = new Matrix(rows);
	for (int i = 0; i < rows; i++) {
		ret->data[i] = data[i*cols + col];
	}
	return ret;
}

Matrix* Matrix::getRow(int row) {
	Matrix* ret = new Matrix(cols);
	memcpy(ret->data, &(data[row*cols]), sizeof(real)*cols);
	return ret;
}

Matrix* Matrix::getRows(int* ids, int n) {
	Matrix* ret = new Matrix(n, cols);
	for (int i = 0; i < n; i++) {
		memcpy(&(ret->data[i*cols]), &(data[ids[i]*cols]), sizeof(real)*cols);
	}
	return ret;
}

Matrix* Matrix::addColumn(int c, real alpha, Matrix* x) {
	if (x->rows != rows) Utils::error("size not match");
	for (int i = 0; i < rows; i++)
		data[i*cols + c] += alpha * x->data[i];
	return this;
}

Matrix* Matrix::addColumn(int c, Matrix *x) {
	return addColumn(c, 1, x);
}

void Matrix::free(Matrix** &As, int n) {
    if (As == NULL) return;
	for (int i = 0; i < n; i++)
		if (As[i] != NULL)
			delete As[i];
	delete[] As;
    As = NULL;
}

void Matrix::free(Matrix*** &As, int n1, int n2) {
    if (As == NULL) return;
	for (int i1 = 0; i1 < n1; i1++) {
		if (As[i1] == NULL) continue;
		for (int i2 = 0; i2 < n2; i2++)
			if (As[i1][i2] != NULL) delete As[i1][i2];
		delete[] As[i1];
	}
	delete[] As;
    As = NULL;
}

Matrix* Matrix::zeros(int rows, int cols) {
	Matrix* ret = new Matrix(rows, cols);
	for (int i = 0; i < ret->length; i++)
		ret->data[i] = 0;
	return ret;
}

Matrix* Matrix::zeros(int dim) {
	return zeros(dim, 1);
}

Matrix* Matrix::uniform(int rows, int cols, real a, real b) {
	Matrix* ret = new Matrix(rows, cols);
	for (int i = 0; i < ret->length; i++) {
		ret->data[i] = rand() * (b-a) / RAND_MAX + a;
	}
	return ret;
}

Matrix* Matrix::uniform(int rows, real a, real b) {
	return uniform(rows, 1, a, b);
}

